Introduction

This document demonstrates how to use the soilDB package to access “soil taxa” mapping (800m grid resolution, CONUS only for now), as part of the ISSR-800 project. Taxa grids were developed from the current SSURGO snapshot using a 800x800 meter grid. Taxa proportions were computed from the geometric intersection of 800m grid and SSURGO polygons (map unit area) and associated component data (component percentage).

Setup

Get required packages.

# run these commands in the R console
install.packages('soilDB', dep = TRUE)
install.packages('mapview', dep = TRUE)
install.packages('rgdal', dep = TRUE)
install.packages('viridis', dep = TRUE)
install.packages('raster', dep = TRUE)
install.packages('rasterVis', dep = TRUE)
install.packages('mapview', dep = TRUE)
install.packages('SoilTaxonomy', dep = TRUE)
# load required libraries
library(soilDB)

library(raster)
library(rasterVis)
library(viridis)

library(rgdal)
library(mapview)

library(sp)
library(spData) 
library(sf)

Prepare state boundaries for CONUS.

# CRS defs used by SEE grids
crs_lower48 <- "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

# use lower '48 state borders as context
# transform to the CRS of the SEE taxa grids
data("us_states")
us_states <- as(us_states, 'Spatial')
us_states <- spTransform(us_states, CRS(crs_lower48))

# get a bbox for CONUS and expand by 100km = 100,000m
b <- bbox(us_states)
x.lim <- c(b[1, 1] - 1e5, b[1, 2] + 1e5)
y.lim <- c(b[2, 1] - 1e5, b[2, 2] + 1e5)

Examples

All vertisols.

taxa <- 'vertisols'
x <- taxaExtent(taxa, level = 'order')
a <- aggregate(x, fact = 5)

## CONUS + states
## cropped to CONUS
levelplot(
  a,
  # maxpixels = 1e5,
  main = names(a),
  margin = FALSE, 
  xlim = x.lim, ylim = y.lim,
  scales = list(draw = FALSE), 
  col.regions = viridis,
  panel = function(...) {
    panel.levelplot(...)
    sp.polygons(us_states, col = 'black', lwd = 2)
  }
)

Xeralfs.

taxa <- 'xeralfs'
x <- taxaExtent(taxa, level = 'suborder')
a <- aggregate(x, fact = 5)

## CONUS + states
## cropped to CONUS
levelplot(
  a,
  # maxpixels = 1e5,
  main = names(a),
  margin = FALSE, 
  xlim = x.lim, ylim = y.lim,
  scales = list(draw = FALSE), 
  col.regions = viridis,
  panel = function(...) {
    panel.levelplot(...)
    sp.polygons(us_states, col = 'black', lwd = 2)
  }
)

Durixeralfs.

taxa <- 'durixeralfs'
x <- taxaExtent(taxa, level = 'greatgroup')
a <- aggregate(x, fact = 5)

## cropped to raster extent
levelplot(
  a,
  # maxpixels = 1e5,
  main = taxa,
  margin = FALSE, 
  # xlim = x.lim, ylim = y.lim,
  scales = list(draw = FALSE), 
  col.regions = viridis,
  panel = function(...) {
    panel.levelplot(...)
    sp.polygons(us_states, col = 'black', lwd = 2)
  }
)

Abruptic durixeralfs.

taxa <- 'abruptic durixeralfs'
x <- taxaExtent(taxa, level = 'subgroup')
a <- aggregate(x, fact = 5)

## cropped to raster extent
levelplot(
  a,
  # maxpixels = 1e5,
  main = taxa,
  margin = FALSE, 
  # xlim = x.lim, ylim = y.lim,
  scales = list(draw = FALSE), 
  col.regions = viridis,
  panel = function(...) {
    panel.levelplot(...)
    sp.polygons(us_states, col = 'black', lwd = 2)
  }
)

Cecil soil series.

s <- seriesExtent('cecil', type = 'raster')
a <- aggregate(s, fact = 5)

## cropped to CONUS
levelplot(
  a,
  # maxpixels = 1e5,
  main = 'Cecil',
  margin = FALSE, 
  # xlim = x.lim, ylim = y.lim,
  scales = list(draw = FALSE), 
  col.regions = viridis,
  panel = function(...) {
    panel.levelplot(...)
    sp.polygons(us_states, col = 'black', lwd = 2)
  }
)

Interactive mapping from R

This is based on the mapview package.

x <- taxaExtent('xeralfs', level = 'suborder')
a <- aggregate(x, fact = 5)

mapview(a, col.regions = viridis, na.color = NA, use.layer.names = TRUE)
x <- taxaExtent('durixeralfs', level = 'greatgroup')
a <- aggregate(x, fact = 5)

mapview(a, col.regions = viridis, na.color = NA, use.layer.names = TRUE)
s <- seriesExtent('san joaquin', type = 'raster')
a <- aggregate(s, fact = 5)

mapview(a, col.regions = viridis, na.color = NA, use.layer.names = TRUE)

Aggregating and mapping multiple soil taxa

This example demonstrates a way of using pattern matching on a soil taxonomic level to filter a list of soil taxa that can be aggregated into one raster and mapped. For example, if we wanted to see all the Andisols and any Andic or Vitrandic subgroups we could use this process. If the pattern match grab too much or isn’t specific enough then you can use the same process to exclude taxa that should not be in the aggregated extent raster (ie. ‘Kandic’ elements).

# load ST data from SoilTaxonomy package
data("ST", package = 'SoilTaxonomy')

# filter to include taxa based on a specified pattern
idx <- grep('and', ST$tax_subgroup)
STgroup <- ST$tax_subgroup[idx]
length(STgroup)
## [1] 573
# filter to remove specified orders and any other based on pattern matching
idx <- grep('ids|ods|ults|ox|els|kandi', STgroup)
STgroup <- STgroup[-idx]
length(STgroup)
## [1] 309
# inspect
head(STgroup)
## [1] "aquandic albaqualfs"     "aquandic endoaqualfs"    "aquandic epiaqualfs"    
## [4] "vitrandic glossocryalfs" "andic glossocryalfs"     "vitrandic haplocryalfs"
# create a list
m <- list()

# loop thru list and fetch taxa extent rasters
for (i in 1:length(STgroup)) {
  m[[i]] <- taxaExtent(x=STgroup[i], level='subgroup', timeout = 60)
}

# merge list of rasters into one rasterlayer
x <- do.call(merge, m)

# give merged raster a name
names(x) <- 'Andisols and Andic Subgroups'

# slippy map
if(requireNamespace("mapview")) {
  mapview::mapview(x, col.regions = viridis::viridis, na.color = NA, use.layer.names = TRUE)
}

This document is based on soilDB version 2.5.9.